arXiv:1506.01689vl [gr-qc] 4Jun2015 


Improvements to the construction of binary black 
hole initial data 

Serguei Ossokine^’^, Francois Foucart^’^, Harald P. Pfeiffer^’^, 
Michael Boyle^, Bela SzilagyP 

^Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, 
Ontario MSS 3H8, Canada 

^Department of Astronomy and Astrophysics, 50 St. George Street, University of 
Toronto, Toronto, ON MSS 3H4, Canada 

^Lawrence Berkeley National Laboratory, 1 Cyclotron Rd, Berkeley, CA 94720, USA; 
Einstein Fellow 

"^Canadian Institute for Advanced Research, 180 Dundas St. West, Toronto, ON MSG 
1Z8, Canada 

^ Center for Radiophysics and Space Research, Cornell University, Ithaca, New York, 
14853 

^Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, CA 
91125 

Abstract. 

Construction of binary black hole initial data is a prerequisite for numerical 
evolutions of binary black holes. This paper reports improvements to the binary black 
hole initial data solver in the Spectral Einstein Code, to allow robust construction of 
initial data for mass-ratio above 10:1, and for dimensionless black hole spins above 
0.9, while improving efficiency for lower mass-ratios and spins. We implement a more 
flexible domain decomposition, adaptive mesh refinement and an updated method for 
choosing free parameters. We also introduce a new method to control and eliminate 
residual linear momentum in initial data for processing systems,and demonstrate that 
it eliminates gravitational mode mixing during the evolution. Finally, the new code is 
applied to construct initial data for hyperbolic scattering and for binaries with very 
small separation. 
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1 . Introduction 

Almost a century ago the existence of gravitational waves was first predicted [T]. 
Gravitational radiation offers an exciting new observational window m E] and the 
enticing possibility of multimessenger astronomy. With the second generation of 
gravitational wave detectors poised to come online [n m E], it is more important than 
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ever to model the likely sources of gravitational waves. Some of the most promising are 
binary black holes, with predicted detection rates of 0.4 — 1000 per year for Advanced 
LIGO [7]. To detect such systems, matched filtering techniques must be used in order to 
extract the signal from the noise [8]. This requires accurate models of binary black hole 
inspiral, merger and ringdown. A vast amount of work has been done in this direction 
in full numerical relativity which is necessary to describe the very dynamic plunge and 
merger regimes (see e.g. PIiniElElDS] for overviews of the field). While many groups 
now successfully simulate binary black hole systems [HI ca El El EH], much of the vast 
7-dimensional parameter space consisting of the mass ratio q and the dimensionless spins 
Xa b remains unexplored. Most of the attention has been focused on binaries close to 
equal mass (g < 8) and modest spin [xa^b 0.8) (although see [El EHl EH |22l |23l l24] i 
For stellar mass black hole binaries, one can expect mass ratios < 15 and arbitrary 
spin magnitudes and orientations, which leads to precession of the spins and the orbital 
plane. Precessing, high mass-ratio binaries have interesting dynamics, causing large 
modulations of the gravitational waveform. One can expect even higher mass ratios 
(g ~ 30) for neutron star-black hole (NSBH) binaries (see [25| a BH-Wolf-Rayet system 
with BH mass 30Mq). At high mass ratios, BBH systems can be used as proxies 
for NSBH systems(e.g. [26]). One would thus like to simulate high-mass ratio BBH 
systems. 

Intermediate mass black holes (IMBH) with masses m = 10^ — IO^Mq have been 
hypothesised to exist to complete the BBH mass hierarchy (e.g., the review [27]b 
Searches for IMBH have been performed and several candidates have been identified 
(see e.g. [281 E9] for recent observations). Higher mass ratio (10 < g < 100) systems 
may serve as models for binaries containing an IMBH and a stellar mass black hole 
or neutron star. Advanced era gravitational wave detectors might be able to observe 
gravitational waves from such systems, with a detection rate of up to 10 events per year 
for stellar-mass - IMBH binaries [7]. It is thus important to explore these systems in 
numerical relativity. 

The first step to numerically evolving a binary black hole spacetime is the 
construction of appropriate data on the initial hypersurface EDI. This involves the 
solution of the elliptic constraint equations with free data that corresponds to a binary 
in quasi-equilibrium, ideally allowing for arbitrary masses, spins and velocities of the 
two black holes. The Spectral Einstein Code (SpEC) [31] includes a BBH initial 
data solver [32] based on the extended conformal thin sandwich equations [33l l34] . 
incorporating quasi-equilibrium black hole boundary conditions [351 [361 ET]. This 
solver has been used to construct BBH for a wide range of configurations [5H] . 
Construction of BBH with increasing mass-ratio, increasing spin magnitudes and the 
desire to construct initial data for highly spinning BBH with arbitrary spin axes have 
necessitated a variety of improvements to the initial data code compared to its original 
presentation [32l [36l |39l HOl SI] . 

This paper summarizes these improvements and extends the original code even 
further, in anticipation of future desire to study even more generic BBH systems. 
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Specifically, here, we present: 

(i) Flexible domain-decomposition to allow a wider range of mass-ratios, spins and 
separations. 

(ii) Adaptive mesh-refinement to enhance computational efficiency and to ensure robust 
numerical convergence for mass-ratios q ^ 5 and dimensionless spins > 0.9. 

(iii) Improved updating formulae for iterative determination of the free parameters. 
These formulas allow one to achieve very high spins and mass ratios, for example an 
equal-mass binary with aligned spins of 0.9999, and a g = 50 single-spinning binary 
with spin of 0.95 on the large black hole. 

(iv) Building on previous work [121113], we control of the ADM linear momentum to 
avoid drifts of the center of mass in BBH evolutions. This eliminates gravitational mode 
mixing due to the motion of the centre of mass with respect to a hxed extraction sphere. 

(v) Control of the center of mass. 

This paper is organized as follows. In Sec. we describe in detail the numerical 
enhancements and additions to the code. In Sec. [^we present the results of initial data 
construction for several challenging conhgurations as well as an exploratory evolution 
of a new data set that demonstrates that the control of linear momentum in initial 
data leads to the elimination of gravitational wave mode mixing. Finally we summarize 
the results in Sec. j^and introduce the construction of initial data for closely separated 
binaries and binaries on hyperbolic orbits as applications of the techniques developed 
in this paper. 

2. Numerical techniques 

The main task of constructing initial data is twofold: hrst, to solve the elliptic constraint 
equations on the initial hypersurface; and then, to ensure that the solution represents the 
astrophysical situation of interest (in our case, a black-hole binary in quasi-equilibrium). 
In SpEC, the former is achieved by using a pseudo-spectral multidomain method; see [32] . 
The number of subdomains is kept hxed, but the resolution of each subdomain is 
dynamically adjusted to obtain low truncation error. To enforce quasi-equilibrium 
conditions, SpEC employs the extended conformal thin sandwich (XCTS) formalism [M] . 
Before solving the conformal thin-sandwich equations, various free parameters must be 
chosen - for example, the sizes of the excision regions, and certain other parameters 
that affect the location, spin or motion of the black holes. The free parameters differ 
from the physical parameters one desires to control, such as the masses and spins of the 
black holes, or the linear momentum Padm of the initial data hypersurface. Therefore, 
iterative root-hnding is needed, as described in Buchman et al jH]. To minimize the 
computational cost associated with many iterations of high resolution solves, we adopt 
a hybrid approach. The resolution of the domain and the free parameters are adjusted 
simulateneously based on the current estimated truncation error and the differences 
between the desired and obtained physical quantities. 

In the remainder of this section, we describe in detail the improvements to the 
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Figure 1. Schematic of the domain decomposition for the initial data solver. The 
thick black circles denote the inner and outer boundaries of the inner spherical shells 
(labeled A and B next to their centers). The blue shaded regions represent five open 
cylinders with axis along the line connecting A and B. The green solidly filled regions 
represent three domains with square cross-section. The thin black circle represents 
the inner boundary of the outer spherical shell, with center indicated by the letter 
C. Dashed lines are guides to the eye, to indicate the dimensions of the various 
subdomains. 


initial data code. 

2.1. Domain decomposition 

Figure [^indicates the geometry of the domain-decomposition employed here. There are 
two inner spherical shells (thick black circles labeled A and B), which are surrounded 
by a set of cylinders (light blue) that are aligned with the axis connecting the two black 
holes. 

Along the axis of the cylinders there are three subdomains with rectangular cross- 
section (indicated in green). One of these is located between the two excision spheres, 
and is a truncated square pyramid. The other two are rectangular blocks. In earlier 
work [32] the two inner spherical shells were restricted to have the same outer radius, 
and all cylinders were restricted to have the same inner radius. This restriction results 
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in a comparatively larger shell around the smaller black hole (B). For very unequal 
mass systems, -C in particular, it may be preferable to have a smaller outer 
radius of shell B, roughly comparable with the sphere of influence of black hole B. This 
would maximize the agreement of the geometry of the domain decomposition with the 
structure of the solution. Therefore, here, we allow unequal radii of the two inner shells, 
as indicated in Fig. [T] This has the largest impact when we consider small separations 
in initial data (for example, for studying remnant properties) where the old domain 
decomposition requires a larger separation between the two black holes than the new 
domain-decomposition in order for the solver to converge. 

The new domain-decomposition uses several parameters from which the placement 
and dimension of each subdomain follow unambiguously. We begin by specifying the 
inner and outer spherical shells: 

• The centres of the inner spherical shells, ca and cb, and of the outer spherical shell, 
Cc. Note that Cc is not required to lie on the line connecting ca and Cb. 

• The inner and outer radii of the inner spherical shells and the outer spherical shell, 
rA, rB, rc, and Ra, Rb, Rc- 

The remaining parameters a, fcyi, /block, and /c determine the relative sizes of the 
cylinders and rectangular blocks: 

• The rectangular blocks and cylinders end on planes orthogonal to the axis 
connecting the centers of the excision spheres. The location of these planes is 
determined by the parameter a, through the requirement that these planes intersect 
the inner spherical shells A and B in circles of radius Ra,b sin a. The opening angle 
of these circles as viewed from the center of the spheres is chosen to have the same 
value for all four planes. 

• The inner radii of the cylinders are determined by the parameter f^yi via 

Pa,b = fcyiRA,B sin a. ( 1 ) 

Note that fcyi < 1 is required for the cylinders 1 and 3 to cover all volume outside 
the spheres A and B. 

• The size of the blocks orthogonal to the line connecting the two spheres is 
determined by the parameter /block, 

(^A,B = fhlockRA,B sin a. ( 2 ) 

The multiplier /block must satisfy /block > fcyi to ensure that the blocks cover the 
entire open region within cylinders 0, 2, and 4. 

• The multipler /c, which measures how much larger the outer size of the cylinders 
is compared to the inner edge of the outer spherical shell: 

ac = Pc = fcr c- ( 3 ) 

To ensure complete overlap between the cylinders and the sphere C, fc > I+C'a/tc, 
with C± being the distance from point C to the axis of the cylinders. 
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The value of /block will determine the relative size of the face of the blocks to the inner 
spheres: If /block > 1, then the edge of the block will be entirely outside the inner 
spherical shell. Conversely, if /block < then the face of the rectangular block is 

completely contained within the inner spherical shell. These considerations will impact 
which subdomain (sphere of cylinder) will provide boundary data for the blocks. 

Our standard values for the grid-internal geometry coefficients are a = 45°, 
/cyi = 0.95, /block = 1-05, and fc = 1.1. We have found these choices to be robust 
for a wide variety of component masses, spins and separations. 

2.2. Adaptive mesh refinement 

An important factor in efficiently generating high-accuracy initial data is the choice of 
resolution in each of the subdomain used in our domain decomposition (see Fig. [^. 
Typically, we want our representation of the solution to have about the same accuracy 
in all subdomains. Unfortunately, we do not know a priori what resolution is needed 
in a given subdomain to reach a target accuracy. Furthermore, the optimal resolution 
varies significantly with the physical parameters of the binary. The old initial data 
solver [321 ES] used hard-coded resolutions, tuned to equal-mass low spin BBH. For 
unequal mass systems, rapidly spinning black holes, and/or widely separated binaries 
the old resolutions are less efficient and can even prevent convergence of the elliptic 
solver when a high accuracy is requested. 

To generate initial data, we generally go through multiple intermediate solves, 
progressively improving the accuracy of the solution while converging towards the 
desired binary parameters. So instead of predetermining the resolution which will be 
used in each subdomain at each level of refinement, we can use the preceeding numerical 
solution to predict the optimal resolution in each subdomain to reach a target accuracy. 
This signihcantly improves the efficiency of the initial data solver, with computing times 
decreased by about an order of magnitude for challenging configurations. And it also 
allows us to push the binary parameters to more extreme values. 

Our multi-domain spectral solver represents the solution in each subdomain as 
a tensor-product of basis-functions. Depending on the topology of the subdomain, 
the basis functions are Chebyshev polynomials, and/or Fourier series, and/or spherical 
harmonics (see [32] for details). 

Following Szilagyi m, for a given subdomain and a given basisfunction, we 
define the power Pi in the i-th mode by the root-mean-square value of all the 
coefficients of the Ath mode across all spectral coefficients of the other basis- 
functions. For instance, in a spherical shell with spectral expansion of the form 

u{r,9,(j)) = Y.!io^Y.i<L,\rn\<i^iirnTiir)Yirn{9,(l)), the radial power would be P* = 
/ _ \ U2 ’ 

(]^ , where = {L + 1)^ represents the number of angular 

coefficient^ 

f For spherical harmonic basis-functions, the top two modes are filtered [32] and are therefore not 
included in the data Pi. 
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For the expected spectral convergence, Pi shonld decay exponentially as a fnnction 
of i im |32], i.e. logj^gP* when plotted vs. i should be a straight line. The slope f of 
this line represents the decrease in the magnitude of the spectral coefficients when going 
from mode i to mode i + 1. We estimate /' using Eq. (53) of Szilagyi [Hj.The current 
truncation error of the spectral expansion is approximated as the highest retained 
coefficient 051 . 

Given the current estimate of the error as e and the estimate of the convergence 
rate as /', we can reach a target accuracy e* by adding 

AN = -^-^ (4) 


modes to the spectral expansion (recall /' < 0 and a higher accuracy means a 
lower e).The answer is rounded up so that AN > 1 if the current accuracy is worse 
than the target accuracy, and we set AN = 0 if e < e*, i.e. the resolution is not allowed 
to decrease. For the conhguration q3 from Table [^the resolution was allowed to decrease 
without noticable impact on the convergence behaviour, cf. Figure 

The outer spherical shell needs comparatively small angular resolution ~ 10, and 
sometimes AMR yields the same resolution at neighboring Pdt- Because the ADM- 
quantities are exclusively evaluated in the outer spherical shell (cf. Sec. 2.4 below), 
this would result in apparent non-convergence of ADM linear and angular momentum. 
Therefore, we increase the angular resolution of the outer sphere by one extra grid- 
point in the 6 direction and the corresponding two extra grid-points in the 0 direction, 
whenever AMR triggers an adjustment to the domain decomposition. 


2.3. Iterative determination of free parameters 

When constructing initial data, we wish to achive desired masses and desired 

black hole spin vectors x*a X*b- The free data, however, is instead given by the radii 
and angular frequencies of the apparent horizons r^^s and which we write as 

u = (rA,rB,rif,rif). (5) 

Therefore, one needs to determine values of the free parameters that result in the desired 
physical parameters. Thus we must solve the system of equations 

F = {Ma - MX, Mb - M*, Xa - Aa, Xb “ X*b) = 0- (6) 

The standard approach to the problem would be to use Newton’s method; however, 
evaluating the Jacobian is too expensive numerically as every evaluation of the 
function F requires an elliptic solve. We instead use the following approach: make 
an initial guess Uq based on the Kerr expressions for both black holes, 

Ma,b = 'i^a,b/{7 + \Jl — 4?’a,s^a,s)> 

Xa,b — ~ ‘^'^a,b^xb^ 


(7) 

( 8 ) 
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and perform an elliptic solve for Fq. We then construct an analytic Jacobian Jq by 
using Eqs. (^) to evaluate the partial derivatives, and update the initial guess by 
= Ho ~ After this we update the Jacobian using Broyden’s method [16] : 

Jk = Jk-i + ^^F(ufc)Au^, (9) 

where Au ^ = This corresponds to the “secant” approximation for a function 

of one variable. Finally we set 


Hfc+i — Hfc Jk 


( 10 ) 


The major advantage of this approach lies in the use of numerical information in 
the update of the Jacobian. This is important in the regime where the simple analytic 
Jacobian becomes inadequate. Broyden’s method is applied to the intrinsic physical 
properties of each black hole, i.e. the eight parameters listed in ([^. We also control more 
general properties of the binary, such as the total linear momentum and the position of 


its centre of mass. As discussed in Sec. |2.5| this is done with explicit updating formulae 
that are applied simulateneously at every step of Broyden’s method. 

We are now faced with two intertwined iterations: AMR to tune grid-sizes to a 
desired truncation error; and root-hnding to adjust free parameters to achieve the desired 
physical masses, spins, etc. When the physical parameters are still far away from the 
desired values, very stringent AMR resolution would waste computing time, so we aim 
to tighten the AMR resolution while simultaneously decreasing root-hnding errors. We 
do so by using an overall truncation error target E£,t for AMR. We start with a large 
value for Edt, corresponding to a small grid-size. As root-hnding residuals decrease, we 
will decrease E^t- We proceed as follows: At iteration A: = 0,1, 2,..., we compute two 
measures of progress in root hnding: First, the residual IZk which quantihes how close 
the physical parameters are to their desired values. IZk is simply the rms error in the 
physical parameters: 


IZk = \ - 



(AM^)2 + (AMb)2 
M2 


+ AAa + Axl + 


adm '^ 


M2 


( 11 ) 


Second, the improvement that indicates how quickly root-hnding converges, dehned 
as 

1/2 


Xk = max 


Qk-zQk-2 


\ Ql-iQl 


k > 3, 


( 12 ) 


where Q' = {AMaAMb, || Ax^II, II AxbII, III’/idmII}. 

We monitor 2 conditions: 


(i) Jk < 

(ii) IZk < ^hEdti 

where Edt is the desired truncation error, and cn and ex are tunable parameters. The 
hrst condition assures that the resolution is increased if the root-hnding convergence 












Improvements to the construction of binary black hole initial data 


9 


becomes “flat” (e.g., due to the inability to measure the masses accurately enough at 
the current resolution). The second condition ensures AMR resolution is sufficiently high 
to ensure the physical parameters can be computed more accurately than the current 
IZk, with e-ji being a safety factor. If either condition is satisfied and we have already 
reached our termination truncation error then the initial data construction is completed. 
Otherwise, we divide Edt by a factor of 10 and continue with the next itertion. For all 
cases we have encountered, the choices 67^ = 10^ and ex = 1-5 have proven to be robust. 


2.4- Calculation of asymptotic quantities 


Accurate knowledge of the total energy, linear momentum and angular momentum of the 
constructed initial data sets aid their characterization. Even more important, accurate 
control of the total linear momentum is essential to avoid a drift of the center of mass 
of the binary during long evolutions, cf. Fig. 

We define the linear and angular momenta on a slice S intersecting spatial inhnity 
on the surface S^o using the Arnowitt-Deser-Misner (ADM) prescription. Our initial 
data satify the asymptotic gauge conditions m 
diij 


o • =C»(r-^), 

dx^ 

7 «A', = 0(r-’). 


(13) 

(14) 


needed to remove ambiguities in the definition of the ADM angular momentum, as well 
as the boundary condition on Soo- The old code |32l EH] directly evaluated 

the resulting surface integrals at infinity nasu, 

PAdm = ^ / (A'“ - Kf’) dS^, (16) 

J SoQ 

= e^Jkx\K^'-K'y^')dSl, (16) 

'X Soo 

using extrapolation in powers of 1/r to infinite radius [IH]. is then found to be 

a combination of 1/r^ terms of and ^ combination of 1/r^ terms. The old 

technique, therefore, is very sensitive to small errors in in the outermost sphere of 
our computational domain (the outer boundary is typically located at rout ~ 10^°M) 
and particularly to the presence of constraint violating modes in that sphere. Typically, 
this leads to large errors in P^^^m resolution, and large errors in even at 

our highest resolution. 

Higher accuracy can be obtained by assuming that the constraints are satisfied on 
our computational domain, and utilizing Gauss’ law to recast the surface integrals on 
S'oo as the sum of a surface integral on a sphere Sq located at a smaller radius and a 
volume integral. Utilizing \h(5oo) = 1, we write 
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Here the normal dSj to Sq points into the interior of Vq (e.g. along +f if it is a coordinate 
sphere) and the factor was inserted to eliminate terms with spatial derivatives of 


from Eq. (21). Using the momentum constraint in the absence of sources, 


V,(A' 




-fOR) = 


dxi 


= 0, (18) 


the volume integral can be simplified to 


p* _ 
-'ADM — 


-hi 


dS, - — / G* dV. 


So 


Svr 


(19) 


'Vo 


Here, 


pij ^ 

G 


K-f 




= f+ f>R 




( 20 ) 
( 21 ) 

where are the connections derived from the conformal metric jij. Note that for 
conformal flatness and maximal slicing, G* = 0 and the volume integral disappears 
(see [50]). 

In practice, for conformally curved initial data. The outer spherical shell extends 
to outer radius 10^°M. Therefore, in the numerical evaluation of the volume integral 


in Eq. (19), the volume element associated with the outermost grid-point becomes very 


large and introduces numerical noise. To avoid this, we roll off the integrand G* beyond 


a certain radius Rr 


G* = 


i.e. we replace G* by G* given by 


r < Rr 



r > Rr. 


( 22 ) 


We choose Rc = 1000 max(tc^, w^), where wa,b are the widths of the Gaussians that 
enforce exponential falloff to conformal flatness (cf. Eqs. 45 and 46 of Lovelace et al 
The ADM angular momentum is also rewritten using Gauss’ law as 


4dm - 


yP^>)dSj-f- f {xG<>-yG^)dV, 

Stt Jvo 


(23) 


with cyclical permutations of {x, y, z) yielding the other components. For maximal 


slicing and conformal flatness in Vq, Eq. (23) simplifies to 

1 

Stt 


4dm = i-t V>(xK« - yK-’) dSr 


(24) 


'S'o 


Because Eq. (23) relies on the cancellation of large volume terms, it can be sensitive 


to errors in . Accordingly, we use Eq. (24) using a surface S'o at sufficiently large 


radius such that in Vq the metric is conformally flat and K = 0. 

To illustrate the importance of the transformations applied to the ADM integrals. 


we consider the convergence test for configuration q50. We evaluate Padm using Eq. (15) 


Figure shows the 


and Eq. (19), and we evaluate Jadm using Eq. (16) and Eq. 
results. 

The calculation of Padm is improved by about one order of magnitude when utilizing 
Gauss’ law, whereas Jadm improves by several orders of magnitude. We point out 
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Figure 2. Accuracy of the calculation of Padm and Jadm for two different methods 
of evaluation. We evaluate Padm and Jadm when truncation error Edt = 10“”“^ is 
reached, and plot differences to the next lower resolution n — 1. Data shown for case 
q50 in Table 


that, in order to achieve any convergence for the old Jadm calculation, we had to 
manually increase the radial resolution in the outer sphere by 1 whenever the domain 
decomposition is adjusted. 

We also compute a new diagnostic, the centre-of-mass CcoM of the initial data 
sets using the formalism developed in Ref. |5T]. In conformal flatness, the expressions 
from [51] reduce to 


CcoM — 


Svr E 


ADM 


lim 

R^oo 


^^ndA, 


(25) 


where n is the outward-pointing unit normal, n = r/r. Equation (25) is numerically 
evaluated by expanding the conformal factor d' in a power-series in 1/r. We read 
off the (angle-dependent) coefficient of the 1/r^ term, and expand this coefficient in 
spherical harmonics. Each individual spherical harmonic term can be integrated against 
n analytically, so that the integral (25) collapses to a linear combination of spherical- 
harmonic coefficients. 
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The quasi-equilibrium conformal thin-sandwich formalism to construct binary black hole 
initial data was developed in a series of papers [Ml ETj [Ml IMl HI]- In this formalism, 
one chooses two excised regions (usually taken to be coordinate spheres) with centres 
ca,b, and solves the extended conformal thin sandwich equations [33ll3l] in the exterior. 
Boundary conditions on the excised regions ensure that they are apparent horizons, and 
control the spin of each black hole. The locations and the sizes of the excised regions 
correlate with the position and masses of the two black holes. Orbital rotation is induced 
by the requirement that certain time-derivatives vanish in a frame rotating with orbital 
velocity fio about the orign. One hnally incorporates a radial expansion factor do, which 
allows hne control of the orbital eccentricity [Ml ISSl [53l SI]- By a suitable choice of 
the conformal quantities, the quasi-equilibrium approach can generate initial data with 
black hole spins of order 0.9998 [M] . 

One shortcoming of the formalism presented in [H] lies in a lack of control of the 
center of mass of the binary, and only incomplete control of the ADM linear momentum 
Padm- The past implementations use the location of the black holes to partially control 
Padm- Consider a small displacement 6c applied to the centres of both excision regions. 
Through the orbital rotation fio about the origin, the displacement Sc induces a change 
in velocity of the black holes of fio x 6c, with a corresponding change in Padm- Therefore, 
Sc could be used to cancel the components of Padm orthogonal to fio; however, the 
cross-product in fig X Sc prevented any correction parallel to fio- For head-on collisions 
with fio = 0 , no control of Padm is possible at all. For the non-precessing simulations 
presented in [H] , the component of Padm parallel to fig vanishes by symmetry, and no 
problems arose. However, for generic precessing binaries, there will be a non-zero linear 
momentum orthogonal to the orbital plane, which results in a drift of the center of mass 
for very long simulations (see [M] for an extreme example). 

Here, we propose a different means to control the full Padm, while simultaneously 
allowing us to control the center of mass as well. We fix the relative separation of the 
centres of the excision spheres, 

ca-cb = D, (26) 


where the separation vector D is user-specified. We use the choice of ca to control the 
center-of-mass CcoM of the binary. Once a first initial data set is computed (with, in 
general, CcoM 7 ^ 0), we can update 


OAjfc+l ^A,k C 


CoM,fc 


MA,fcAMB,fc — MB,fcAMA,fc 


D. 


(27) 


(Mi,i + 

With the black-hole centres now used to control the centre of mass, we need a 
different means to control Padm- We add in the outer boundary condition on the shift 
(Eq. (38c) of [iQ]) a constant velocity vq: 


/?* = (f2o X r)* -h doC + Vq 


on B. 


(28) 
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Figure 3. Control of the centre of mass and linear momentum for a g = 10 generically 
processing binary (see case qlO in Table . Shown are the magnitude of CcoM (black 
solid lines) and Padm (red dashed lines) as a function of root-finding iteration k. 


Here B represents the outer boundary, a sphere with radius R = The velocity 

vq will effect the overall motion of the binary, and will be reflected in a corresponding 
change in Padm by -EadmVq, where Eadm is the ADM-energy of the binary. During 
iterative root-finding of the free parameters, we adjust vq to achieve Padm = C[§} 

To motivate the updating formula for vq, consider a perturbation of vq by (5vo, and 
a perturbation of ca by 5c. If we allow the masses to vary, then a Newtonian-inspired 
formula is 


Vo,fc-|-l — Vo,fe- 


ADM,fc 


f(AMA,fc+AMB,fc)(Vfc-briXCA,fc)-rix5cA,fc- 


AMi 


B,k 


Mk ’ ’ ’ ’ Mk 

To summarize, relative to earlier initial-data sets, we modify the outer boundary 
condition for the shift by the term vq, cf. Eq. (28), and use updating formulae (27) 
and (29) to adjust ca and vq. Section 2.4 describes how we compute Padm and CcoM- 
We demonstrate the efficiency of the updating formulas Eqs. ( 27p9 ) in Fig. |^that 
shows the magnitude of CcoM and Padm as a function of root-finding iteration for a 
g = 10 processing binary (case qlO in Table [^. The convergence is evidently very fast, 
with the final values of ~ 10“® and ~ 10“® respectively. This means that even for an 
inspiral lasting 10® M, the drift of the centre of mass due to residual linear momentum 
in initial data will be only ~ O.OIM. 


§ Using the obtained vector /3* as the shift-vector in an evolution results in a translating outer boundary; 
this effect is eliminated by evolving with a shift vector of /3* — Ug. 


rixD.(29) 
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Table 1. Physical parameters of the initial data sets used for testing the new initial 
data code. 


Name 

q 

Xi 

X2 

Do/M 

MQo 

SpinO.9999 

1 

(0, 0, 0.9999) 

(0, 0, 0.9999) 

14.17 

0.01682 

q3 

3 

(0, 0.49, -0.755) 

(0, 0, 0) 

15.48 

0.01515 

qlO 

10 

(0.815, -0.203, 0.525) 

(-0.087, 0.619, 0.647) 

15.09 

0.01547 

qSO 

50 

(-0.045, 0.646, -0.695) 

(0, 0, 0) 

16 

0.01428 


3. Numerical results 


3.1. Initial data construction 

We test the improvements described in the previous sections on several cases of 
interest, whose parameters are summarized in Table The parameters were chosen 
to demonstrate the range of initial data sets that can be constructed with the new code 
and to provide some overlap with regions of parameter space which could be achieved 
previously. 


We hrst illustrate the performance of the AMR outlined in Sec. 2^ with the case 
q3, a conhguration we will compare with the old BBH solver below. To demonstrate 
AMR in isolation, we £x initial data parameters, and start with target truncation error 
Ejjt = 10“^. We solve the constraint equations, estimate spectral truncation errors 
and update numerical resolution via Eq. Q. Whenever we reach the desired truncation 
error, we tighten the AMR error tolerances by dividing by 10, until a truncation 
error of 10“® is reached. Figure [^illustrates the behaviour of the AMR algorithm during 
this test. The top panel shows the total number of collocation points in the domain, 
which grows with each AMR iteration. The bottom panel demonstrates that the largest 
truncation error across all subdomains, maxe, closely tracks the truncation error target 
Edt- 

Figure [^ shows a convergence test of the AMR sequence shown in Fig. [^ Plotted 
are various quantities as a function of the effective number of grid-points The top 

panel demonstrates the exponential decrease in the norms of the Hamiltonian and 
momentum constraints, which implies that this data set is constraint-satisfying. The 
constraints are given explicitly by: 




Ham 






KabK^’^) , 


CMom — DhK^ — DaK, 


where D is the covariant derivative associated with the spatial metric. The 
simply the normalized pointwise norm over all collocation points: 


(30) 

(31) 
norm is 


P = 




N 


N 




(32) 
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n 


Figure 4. Behaviour of the AMR algorithm for case q3 as function of AMR iteration 
n. Top: total number of collocation points. Bottom: Highest truncation error in all 
subdomains, and target truncation error. The free parameters are fixed to their values 
at the end of root-finding. 


The convergence of the masses and spins is shown in the middle panel. Here we 
plot the norms of the differences between the qnantity at a given iteration and its value 
at the highest resolution: 

AQ = IQi - QmaxI- (33) 

Once again, the convergence is essentially exponential. The bottom panel of Fig. 
hnally examines the convergence of the ADM quantities and the center-of-mass 
computation. Though convergence is not as clean as for the constraints, the bottom 
panel of Figure shows that all the asymptotic quantities can be determined to better 
than 10“®. 

To conclude our detailed examination of the initial data set q3, we contrast the new 
code described here with the old code [32l Hi]. One of the most important upgrades 
lies in the control of Padm- Figure shows the components of Padm as a function 
of root-hnding iteration k for both the new and the old code|j|] Both codes successfully 
drive Tadm -^adm expected. But only the new code also drives Dadm 

zero, whereas the old code yields PXdm — 0.00138 . As we shall see in the next section, 
this produces non-trivial differences in the evolution. 

Turning to the more challenging cases listed in Table we have performed similar 
tests to case q3, with the free parameters hxed to their values at the end of root 
hnding and only the resolution changing from iteration to iteration. As an example, 

II Both codes compute Padm in the same way (via ??), but differ in the way it is controlled. 
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Figure 5. Convergence of physical quantities with increasing resolution for case 
q3. Top: the norm of Hamiltonian and momentum constraints, middle: masses and 
spins, bottom: ADM quantities and center of mass. “A” indicates the difference 
between the value of the respective quantity at the current resolution and the highest 
resolutions. The free parameters are fixed to their values at the end of root finding. 


Figure shows a subset of the convergence data. This hgure demonstrates that the 
exponential convergence shown previously for case q3 extends to all cases. In particular, 
the constraints are exponentially convergent. All four cases complete with a maximal 
resolution of less than 2.5 x 10® points, an improvement of a factor of 2 — 3 over the old 
codejf] Owing to the more challenging configurations, however, the constraints are 1-2 
orders of magnitude larger. The physical parameters are also exponentially converegent 
with resolution, as illustrated in the lower panel of Fig. We use AMb since it is 
frequently harder in a high mass ratio to resolve the smaller black hole, so this provides 
a conservative convergence test. 

3.2. Root-finding 

It is also important to examine the performance of the updated root-hnding procedure 
based on Broyden’s method. Figure shows the root-finding results for cases q3, qlO 

^ We note that the case qSO could not be constructed with the old code. 








AM Hamiltonian constraint 
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Figure 6. Comparison of Padm control between the new code (solid lines) and 
the old code (dashed lines) as a function of root-hnding iteration k. Both versions 
successfully control the x and y components of the linear momentum, but only the 
new code controls the z-component as well. 



Figure 7. Overview of initial data results for cases in Table B Top: convergence 
of the norm of the Hamiltonian constraint. Bottom: convergence of the mass 
of the smaller black hole.The free parameters are fixed to their values at the end of 
root-finding. 
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Figure 8. Convergence of the root-finding procedure for masses and spins when the 
old updating formulae are used (dashed lines), and with the new updating formulae 
developed here (solid lines). For q3 both algorithms perform well, whereas for qlO the 
new code converges about twice as fast as the old code. Finally for SpinO.9999 the old 
code fails to achieve desired masses and spins, while the new code gives errors of order 
10 -®. 

and SpinO.9999 done with the old and new versions of the code. Note that dnring root- 
hnding the resolntion of the subdomains is also allowed to change to achieve the desired 
truncation error. For low mass ratio both codes show similar rates of convergence and 
hnal errors. The situation changes for case qlO, where the old code has trouble achieving 
low errors in masses and spins, while the new root-hnding procedure described in Sec. m 
results in errors of order 10 Finally, for case SpinO.9999, the results are drastically 
different: the old code has errors in the masses of order several x 10“^ and spins of order 
10“^. Since we are attempting to construct a binary with dimensionless spins of 0.9999 
it becomes clear that the old code is inadequate for this purpose. On the other hand, 
the new root-hnding procedure successfully reduces the errors in physical quantities to 
the level of 10“®. Thus, the new algorithm allows us to achieve the desired values of the 
physical quantities which is especially important as we push to higher spin magnitudes. 

On average, the new code is about 25-50% as fast as the old one. For example, for 
the case qlO, the old code took 12.4 hours to complete, whereas the new took 6 hours 
on 12 cores of a Westmere node of the Briaree compute cluster. Therefore, the new code 
is indeed more efficient than the old while achieving the same or better accuracy. 

3.3. Exploratory evolution 

We have emphasized above the importance of controlling Padm- We now evolve initial 
data for case q3 constructed with the old and the new initial-data code, and compare 
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Figure 9. Top: Evolution of the normalized constraints. Middle: Evolution of the 
Christodolou mass of the large black hole. Bottom: Evolution of the spin \xa\- 


the two evolutions in detail. 

We being by considering the convergence of constraints and quasi-local quantities 
during the evolution. The top panel of Figure shows the norm of the normalized 
constraint violations during the evolution (see Eq.(71) of [55]). It is obvious that 
both codes show similar convergence properties, as expected. Further, the initial 
spike of constraint violations due to junk radiation is virtually indistinguishable, which 
indicates that the new method of constructing initial data does not introduce additional 
constraint-violating modes. The middle and bottom panels of Fig. [^ show the evolution 
of the Christodoulou mass and the spin magnitude of the large black hole. The 
differences between the evolutions of the old and new initial data sets are consistent 
with truncation error. Thus we conclude that the quasi-local quantities are the same in 
both data sets. 

Turning attention to the trajectories of the black holes, we hnd a stark difference in 
the evolutions. Figure [TO] shows the motion of the large black hole in inertial coordinates 
for both runs. The uncontrolled residual linear momentum Fadm the old initial data 
causes the centre of mass of the binary to drift linearly during the evolution, as shown in 
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Figure 10. Left: The trajectories of the centres of the apparent horizons of the big 
black holes in the intertial frame. The thick black line represents the motion of the 
coordinate centre of mass for the drifting binary. The thin black curve correspond to 
initial data with large drift, the thick red, without. Right: The components of the 
coordinate centre of mass for both runs. The dashed curves refer to initial data with 
large drift, the solid curves without. 


the right panel of Figure W Such a drift may have multiple undesirable consequences. 
Most immediately, it causes the gravitational wave extraction spheres to be off-center 
from the center-of-mass of the binary, which will cause mixing of the spherical harmonic 
modes of the gravitational radiation, an effect discussed in more detail below. Moreover, 
SpEC’s constraint preserving outer boundary conditions [SSI EZl ES] are designed to 
work best for low-order spherical harmonic modes of the outgoing radiation. If the 
binary is offset relative to the outer boundary (for instance due to a drift of the center 
of mass), higher order spherical harmonic components will become more important, 
possibly leading to an additional runaway acceleration of the center of mass [59] . 

To examine the dynamics of the binary, we study the orbital frequency vector 
n = r X r/|rp. The left panel of Figure O shows the projection of Q onto the unit 
sphere, making it apparent that the precession and nutation dynamics are very similar 
until very close to merger. The right panel shows a plot of = |f2| from which several 
features are apparent. The evolution of is qualitatively the same in both cases, 
consistent with expectation that removing a coordinate motion of the centre of mass does 
not change the binary dynamics. Likewise the initial pulse of junk radiation (inset A) 
appears quite similar. However, small oscillations in H are more pronounced in the new 
code (inset C). This is reflected in the measured values of the eccentricities: e = 10“^ for 
the old, e = 2.5 x 10“^ for the new code. The difference in eccentricity arises because the 
new term n* in the outer boundary condition Eq. (28) does slightly modify the content 
of the initial data. In this particular case, |no| ~ 10“^, so that it is not unreasonable to 
expect the orbital eccentricity to change by a comparable magnitude. The initial orbital 
frequency Hq, initial radial velocity do, and initial separation Dq listed in Table [T] were 
tuned to result in essentially vanishing eccentricity in the old initial data [33] • The new 
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Figure 11. Left: The direction 0(t) of the angular velocity vector on the unit sphere. 
Note the excellent agreement in precession dynamics. Right: The magnitude Il{t) of 
the angular velocity vector. The overall agreement is very good. Insets A and B 
highlight the different times to merger due to junk radiation dynamics, while inset C 
demonstrates the different eccentricities. 


initial data constructed from the identical initial data parameters must therefore have 
a slightly larger eccentricity. If we had tuned to vanishing eccentricity with the new 
initial data, then the old initial data would exhibit the larger eccentricity. 

The evolutions of the old and new initial data also result in a different time to 
merger, cf. panel B of Fig. 11 This difference could be caused either by the slightly 


different inspiral dynamics like eccentricity, or could simply be due to truncation error 
of our low resolution evolution. 

Finally, we examine the waveforms for the two runs. Most strikingly, the movement 
of the coordinate centre of mass shown in Figure 10 is also reflected in the spherical- 
harmonic decomposition of the waveform. This is most easily seen in the sub-dominant 
modes. Figure 12 shows the {i,m) = (3,1) modes of the spin-weighted spherical- 
harmonic (SWSH) decompositions of the waveforms measured from the old initial 
data and measured from the new initial data. Since gravitational waves in SpEC 
are extracted on a coordinate sphere centered on the origin, a drifting source mixes the 
modes of the SWSH decomposition. As seen in the lower panel of the hgure, this mixing 
introduces very large effects. To verify that these effects are primarily due to the motion 
of the center of mass, we have also transformed h"®" to a frame in which the center of 
mass is moving as in the original initial data. The initial position of h"®" is transformed 
to agree with the center of mass of the old initial data as measured by Eq. (25), and its 
velocity is transformed to agree with Padm/Mabm of the old initial data as measured 
by Eq. (19). This transformation is applied entirely at future null inhnity by the method 
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[t - r^)/M 


Figure 12. The /13 i waveform modes, as measured in evolutions of the original 
and the new initial data, and extrapolated to future null infinity, . The upper 

panel shows the waveform h"®’" from the new initial data, measured in a frame that is 
centered at the origin of the simulation coordinates. The lower panel shows the same 
data with a transformation applied on as described in the text, as well as the 
waveform from the original initial data measured in its simulation coordinates— 
in which the black holes are moving as shown in Fig. Essentially, the center of 
mass is stationary at the origin in the upper panel, and is moving in the lower panel. 


described in [60], and is a special case of a BMS transformation [6I1162]. It will thus be 
seen in any waveforms, whether extrapolated jM] (as seen here) or extracted by Cauchy 


characteristic methods [Ml ESllMl ET]. As shown in the lower panel of Fig. the 
transformation reproduces the features seen in very well. 

Mode decompositions like this one are used very frequently for analyzing numerical 
models, and for constructing analytical models. If they are unmodeled and uncontrolled, 
effects like those seen in the lower panel will simply appear to be errors in the waveform. 
This could negatively impact uncertainty estimates of numerical simulations, error 
estimates for analytical waveforms, or calibration of waveform models to numerical 
results. These effects will also be present in any calculation that uses the waveforms 
to compute physical quantities such as the flux of linear and angular momentum. By 
removing extraneous displacements and boosts, this new initial data code simplifies such 
analyses]^ 


The drift described here is a linear motion due to residual linear momentum in initial data. 
Controlling this drift will not help for other types of motion present in very long simulations; see |59j . 
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4. Discussion 

Numerical evolution of binary black hole spacetimes requires accurate initial data. In 
this work we have improved the initial data techniques in SpEC to allow access to a much 
wider parameter space of generically processing high mass ratio, high-spin binaries. 
A more flexible domain decomposition allows for stable solution for high-mass ratio 
and high spin binaries. An enhanced root-hnding algorithm is used to achieve desired 
physical parameters for the binary. This becomes important when a naive analytic 
Jacobian is not appropriate, which is precisely the case for high mass ratios and spins, 
see Figure]^ Adaptive mesh rehnement drastically improves efficiency and robustness of 
the code, displaying exponential convergence of constraints, c.f. Figure Finally, a new 
method to control the linear momentum is used to eliminate a linear drift of the centre 
of mass during evolution. This in turn nullihes spurious gravitational mode mixing, 
which is of paramount importance for construction of hybrid waveforms or calibration 
of phenomenological models as demonstrated by Figure [T^ 

An interesting application of the improved initial data code is the construction 
of initial data for hyperbolic encounters. Such systems have been studied in the past 
(e.g. [68116^ 170]) and provide a laboratory for exploring strong held physics in a different 
regime than the binary inspiral. Using the new code, we have successfully constructed 
initial data for hyperbolic encounters for a selection of mass ratios and spins, which was 
not possible before in SpEC. As a simple example, we evolve two systems of two equal 
mass black holes that are initially separated by 60M and have a velocity of ~ 0.14c. 
Both systems have the same impact parameter = 15M, and differ only in the black 
hole spins: In one case the black holes are non-spinning, in the other both holes have 
dimensionless spins x = 0.5 initially in the x direction. Figure [T^ shows the trajectories 
of the two black holes. In the presence of spin, the spin-orbit interactions cause the 
plane of scattering to change and also change the deflection angle of the hyperbolic 
encounter. Exploration of other parameters is left to future investigations. 

Another application is the construction of initial data for binaries with very small 
initial separation, corresponding to only a few orbits before merger. This is useful if 
one is interested in the properties of the merger remnant, e.g. for calibrating analytical 
waveform models but evolving a long inspiral is too computationally expensive. As an 
example, we construct initial data for a system with q = 21, xi = 0.66, X 2 = 0.41 
(oriented in random directions) and initial orbital frequency of MQ = 0.032, and initial 
coordinate separation Do = 8.82 M. We note that initial data for binaries near ISCO at 
high mass-ratio is challenging and further work remains to be done to make it robust 
for g > 10 regime. 
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Figure 13. Hyperbolic encounter of two equal mass black holes shown through the 
coordinate trajectories of the black holes. Left: non-spinning black holes. Right: 
black holes with spin Xi = X 2 = (0.5, 0.0,0.0). The black holes start on the cc-axis in 
the X — y plane (shown in grey). In the spinning case the motion is not confined to 
this plane. 
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